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Abstract. We solve the Bogoliubov-de Gennes equations for an inhomogeneous condensate 
in the vicinity of a linear turning point. A stable integration scheme is developed using a 
transformation into an adiabatic basis. We identify boundary modes trapped in a potential 
whose shape is similar to a Hartree-Fock mean-field treatment. These modes are non- 
resonantly excited when bulk modes reflect at the turning point and contribute significantly 
to the spectrum of local density fluctuations. 
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Introduction 


he achievement of Bose-E instein condensation in ultra-cold trapped atomic gases 


(iPitaevskii and Stringarii 120031) has provided experimentalists with a ‘direct look’ at quantum 


mechanical wave functions. In addition, the atom-atom interactions that become relevant 
despite the low densities, lead to a nonlinear wave mechanics of degenerate B ose gases, as 


described by the celebrated Gross-Pitaevskii equation at the mean field level (IGross . 


1961 


Pitaevskiii Il96lh . see Eq.© below. Nonlinearity brings in qualitatively new features in 
inhomogeneous systems, for example: by neglecting the kinetic energy (second derivative), 
one gets a nontrivial solution with a fixed amplitude, the so-called Thomas-Fermi condensate. 
This approximation breaks down in the vicin ity of a turning point, and the condensate’ s 


kinetic energy acquires logarithmic corrections (IDalfovo et al. . 


1996; 


Fetter and Feder. 


1998 ). 


One has to deal with a nonlinear boundary la yer problem, similar to the G inzburg-Fandau 
description of the surface of a superconductor iFifshitz and Pitaevskii . 1980 ) that leads to the 
distinction between type I and II superconductors. 

^ henkel@uni-potsdam.de 




























Excitations at the border of a condensate 


2 


We address in this paper the wave meehanies of elementary exeitations around the Gross- 
Pitaevskii equation by foeusing on a typieal turning point where the trapping potential is 
approximately linear. This situation is of eourse well known for the linear Sehrodinger 
equation: it leads to an Airy funetion and the famous 7r/4 phase when semielassieal wave 


iLangei, 

1937; 

Messiah. 

1995) 


Messiahl 119951) . In the nonlinear ease, one is dealing with two eoupled wave 


funetions or Bogoliubov-de Gennes (BdG) modes u and v. This eoniplieat es the semielassieal 
analysis and has led to modified WKB teehniques (iHvouguehi et al.l.l2002h . A straightforward 
numerieal approaeh is impossible beeause the higher (fourth) order of the wave equation 
aetually generates an instability. One of the motivations of the present analysis is to provide 
a robust seheme for the BdG modes that ean be used as a stepping stone for inhomogeneous 
low-dimensional Bose gases at finite temperature. Indeed, in this ease, thermally exeited 
modes give a domin ant eontributiqn in the infrared and enforee the introduetion of the quasi- 


eondensate eoneept (iKagan et al.l. l2000l: 


Andersen et al. 


2 OO 2 I: [Mora and CastinL 120031) . The 


Bogoliubov modes that we derive here eapture the role of spatial eoherenee (deloealised 
waves) and may provide a quantitative assessment of the physies beyond the loeal density 
approximation. Indeed, we find that spatial gradients of the eondensate density play a key 
role for the elementary exeitations in the border region. 

The paper is organised as follows. We reeall the mean-field theory for the elementary 
exeitations of an inhomogeneous degenerate Bose gas and formulate the boundary eonditions 
on both sides of the position where the ehemieal potential erosses a linear(ised) trapping 
potential (Seed])- In SeelH the BdG equations are solved approximately with the help of an 
adiabatie basis that generalises the transformation to density and phase modes u ± v. We 
diseuss in partieular the appearanee of ‘trapped modes’ near the eondensate boundary. The 
eonsequenees for physieal observables like the eondensate depletion, t he average thermal 
density and its fiuetuations are illustrated in SeeS In a eompanion paper dpiallo and Henkel 
20151), we analyse the eorreetion to the WKB (Ganger) phase at the nonlinear turning point 
and its role for the speetral density of elementary modes. 


1. Model 


Interaeting Bose gases at low temperatures are quite sueeessfully deseribed by a mean-field 
theory provided most of the partieles oeeupy the eondensate mode. This mode then solves a 
non-linear Sehrodinger equation, also known as the Gross-Pitaevskii equation (GPE): 

- 77— + Vf H - 

2m m 


m 


( 1 ) 
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This is the stationary form of the GPE, with the eigenvalue /x called the chemical potential. 
The (positive) scattering length specifies the density-dependence of the inter-particle 
interactions at the mean-field level, and V = E(r) is an external potential. In this paper, 
we focus on a quasi-one-dimensional trap and replace the interaction term by an effective 
interaction strength proportional to Og and the transverse confinement. In addition, we 
focus on the spatial region where the potential can be linearized, more specifically in the 
vicinity of a turning point: V(z) ~ — Fz. By shifting the ^-coordinate, the chemical 

potential drops out of the G PE. With the proper choice of units (see Table [U), the GPE 
finally takes a universal form iDalfovo et ah . 1996 ). also recognisable as the second Painleve 


transcendent (lAblowitz and 


Seguj, 


19771) . 




Zf + |0f0 = 0 


( 2 ) 


The linearisation around the mean field leads to the Bogoliubov-de Gennes equations that in 
the same units can be written as 


dz^ 

d^v 


dz^ 


ZU -\r 2|0fM -I- (p^V* = Eu 
— zv + 2\(j)\‘^v + cjfu* = — Ev 


(3) 


where E > 0 is the frequency (energy) of the elementary excitation, measured relative to the 
chemical potential. We fix the phases of 0 , m , v to be real, choosing positive 0 . 


System 

Eength 

Temperature 

Erequency 

Density* 


^ = 0^/{2mFY/^ FljkB 

Fi/fi-Kb) 

Ftig 

ID, gravity 

0.3 /im 

31 nK 

640 Hz 

6.4//im 

ID, 100 /im length, 10 Hz trap 

1.1 /im 

2.3 nK 

48 Hz 

0.47//im 

3D, 100 /im diam, 30 Hz trap 

0.53 /im 

9.9 nK 

210 Hz 

28//im^ 


Table 1. Natural units for the Bogoliubov-de Gennes (BdG) equations in a linear potential 
V = —Fz. The interaction constant in the quasi-lD geometries (first and second line) is 
g = 2Tiuj±as with transvers e trapping fre q uency u:±I2tt = 10 kHz and s-wave scattering 


length Us = 95 ag for Rb87 (lEgorov et al 


201 11) . The trapped systems are considered in a 


harmonic confinement, the potential being linearised at the Th omas-Fermi radius. 


*The density scale for the 3D trap is taken as l/(87rasf^) (IDalfovo et al 


1996). 


1.1. Condensate wave function 

The physically relevant solution to Eq.(l2l) is known as the second Painleve transcendent and 
interpolates smoothly from an exponentially decreasing (tunnelling) wave to the Thomas- 
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Figure 1. Condensate wave function (second Painleve transcendent and solution to Eq.©) 
and its asymptotic behaviour (Eqs.(|4l|5]l, dashed). We keep only the leading term in Eq.(|5]l. 


Fermi solution obtained by neglecting the second derivative (Fig[T]). Since one deals 


arbitrary, and it has 

)een shown that ( 

Ablowitz and Segurl. 

1977; 

Hastings and McLeod 

1980; 

Dalfovo et al. 

1996 

:Lundh et ah. 

1997 

) 


2; —>■ —cxo 

On the dense side, 


(j){z) —)■ \/2Ki{—z) 


(4) 


Lundh et al. 


(1 19971) and lMargetisI (I 2 OOOI) have improved the Thomas-Fermi 


solution into the expansion 

2; —)■ +CXD : 


(t){z) V^(l 


Cl 


C2 


(5) 


^hf(^) — —z + 2|0p — 


z3 ^6 ) 

with coefficients Ci = 1/8, C 2 = 73/128, ...The condensate density appears in the BdG 
Eqs.dl]) for u and v, for example via the Hartree-Fock potential 

-z forzC-l 

(b) 

+z for z +1 

which is a wedge-shaped trap (Fig|2t^c/0, thin solid line). Numerically, we find a polynomial 
approximation near its bottom (error < 0.01) 

Vhf(^) ^Vo + V2(z - Zof + Vsiz - Zo)^ (7) 

for jz — ^o| < 1 


with the minimum located at uq ~ 0.53 and zq ~ 0.13 and parameters V 2 ~ 0.47, V 3 ~ 0.041. 
It will turn out, however, that the Hartree-Fock well is irrelevant for the Bogoliubov solutions 
- the only message to keep is the characteristic energy scale E = (9(1). 


1.2. Boundary conditions for Bogoliubov solutions 

The Bogoliubov modes feature an intermediate zone —E<z<E where the excitation 
changes its character from ‘single-particle’ to ‘collective’. Outside this zone, the asymptotic 
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behaviour is as follows. 

On the dilute side, the eondensate 0(z) in Eqs.® vanishes, and the mode functions u 
and V decouple. The linear branch of the Hartree-Fock potential V}i-p{z) — 2 ; is a good 
approximation. We thus have a turning point ze for u{z) at ze ~ —E. The mode v{z) 
is already in the tunneling regime for ^ < -1 because of the opposite sign of the energy 
eigenvalue in Eq.®. 




Figure 2. (left) Illustration of the ‘Hartree-Fock potential’ (Eq.®, thin solid line, middle) and 
its variants that appear in the equations ® for ‘phase’ (blue, bottom) and ‘density’ (red, top) 
modes. Dashed; parabolic approximation (|7]l to the Hartree-Fock potential. 

(right) Potentials in the adiabatic approximation for three energies. Upper (red) curves: 
‘density mode’ k^{z), lower (blue) curves: ‘phase mode’ —k'^{z) (see Eas. dT^ fT^ k As 
the energy E or the condensate mean-field potential (f{z) increases, the potentials are pushed 
apart. The physical mode functions correspond, in this representation, to solutions at zero 
energy (thick black line). The bump around — 1 < z < 0 at low energies is due to the 
geometric potential, see discussion after Eas. lfTSlfT^ below. Upper dashed line; Hartree-Fock 
potential (Eq.®, see left panel), lower dashed line: Coulomb-like asymptote of Eq.® . 

In the dense region where the condensate dominates, also the coupling ~ fizY between 
u and V becomes large. It is convenient to switch to the ‘density-phase’ representation 
/ = (u + v)/ \/2 and g = {u- v) /\/2 where the equations become 

- ^)9 = Ef ( 8 ) 

The ‘potentials’ that appear here are plotted in red and blue in Fig®Ze/t). The ‘density mode’ 
/ corresponds to a well (upper red) whose spectrum starts above zero energy (the minimum 
value of 30^ — is ~ 0.78 at z ^ —0.21). It is ‘enslaved’ to the ‘phase mode’ g that appears 
as a source term Eq.([8l), first line. The potential for the phase mode is a smooth barrier 
that crosses zero at z ~ 0.8 and vanishes for z S> 1 (Fig® left), lower blue). To take into 
account the density-phase coupling proportional to E, we perform the adiabatic elimination 
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/ ~ Eg/{2z), using the Thomas-Fermi asymptote 30^ 


2z and neglecting the second 


derivative. This gives deep in the condensate the equation for the phase mode 


2 ; —)■ CX3 


d‘^g ^2 


dz^ 2z^ 


0 


(9) 


This one-dimensional Coulomb problem has exact solutions that are discussed in Sec l2.2l 
below. To state the boundary conditions, a simpler semiclassical (WKB) treatment will 
suffice. From Eq.®, identify the local wavenumber k{z) = E/sf^ and calculate the action 
integral: one gets two independent solutions from the real and imaginary parts of 

( 2 zfl^ 


2 ; —)■ cxo : 






■ exp(iE\/^) 


( 10 ) 


Since f{z) is smaller by a factor E/{2z), this expression will dominate the behaviour of 
both u{z) and v{z) deep in the condensate. We call this asymptote the ‘local density 
approximation’ because the WKB treatment assumes that the condensate density (/“^{z) ~ z 
varies slowly enough. In terms of the wave number, we require \dk/dz\ -C or Ey/^ ^ 1. 
This condition illustrates that the border region 2 ; ~ 0 and the low-energy range E -C 1 are 
actually challen ging and require technique s beyond the WKB approximation. For a discussion 
of this point, see lPiallo and Henkell (120151) . 

To summarise, the physically relevant boundary conditions are 


^ <C -1 : 


( 11 ) 


(i) dilute domain through the turning point z ~ —E, but away from the condensate border 

u{z) = aPd[—E — z) 
v{z) = pA.i{E-z) 

This covers the tunnelling region where both Airy functions become exponentially small. 
At large energies, the solutions are such that v{z) is much smaller than u{z). 

(ii) local density approximation in the dense (condensate) region 

z:^ E, 1/E^ : 

u{z) = ^ cos{Es/^ — 7r/4 -f 5) 

(2^)V4 


v{z) = 




cos{Ey/^ — 7r/4 -f 5) 


( 12 ) 


We have considered here only the leading order terms proportional to the phase mode 
g (Eq. tfT^ l. The normalisation is such that the solutions 5 = 0,7r/2 have the same 
amplitude and unit Wronskian, see Appendix A for details. 

The phase shift 5 in the bulk asymptote (fT^ depends on the relative weight between real and 
imaginary parts of the complex solutions (fTOl) . The reference —7r/4 is explained in Sec J2.2[ 
We emphasise that this phase shift 5 = 6{E) ‘carries’ information about the behaviour near 
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the condensate border into the bulk. For a matching of the BdG solutions near the turning 
point w ith bulk solutions using boundary layer techniques, see for example 
il998L 


Fetter and Fedei 


2. Phase and density modes 

The coupled BdG equations contain unphysical solutions that grow for 2 ; ±00 and 

that typically contaminate numerical trials when the BdG equations are straightforwardly 
integrated. This can be seen from the second line of Eqs.dS]) whose homogeneous solutions 
are ‘under the barrier’ and grow exponentially. We have developed instead a semi-analytical 
scheme where the unstable modes are eliminated. The idea is to perform a rotation in the 
MU-plane that diagonalises the coupling. 


2.1. Adiabatic transformation 


We make the following Ansatz for a rotated set of amplitudes 


/ M \ / COS0/2 sin6*/2 \ \ 

\ u / \ —sin0/2 cos 0/2 / \ u / 


(13) 


and find that the coupling between u and v (Eqs.dS])) is removed when the rotation angle 0 is 
chosen as 

tan0(z) = ^^ (14) 

Note that in the dense region, we have 0 —)■ 7r/2 and the amplitudes u, v approach the phase 
and density modes {g, f) introduced above Eq.dS]). We note that a hyperbolic rotation that 
preserves the Bogoliubov norm can also be used, but leads only to minor changes in 

notation. 

The equations for u and v do not decouple completely because the rotation angle 0 
depends on position. By working out the second derivative of Eq.([T3l). we get 


(Pu 
dz'^ 
d^v 


= Lv 


dz'^ 


+ = — Lu 


(15) 


(16) 


where the coupling involves derivatives of the condensate density via the differential operator 


L = \e" + 0'^ (17) 

The ‘adiabatic potentials’ —P and in Eqs. dTSlfT^ are recognised as the generalisations of 
the phase and density potentials (jf — z and 30^ — zof Eqs.®. They are plotted in Fig^right) 
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-e= -z + 2(f^-^E^ + cl)^+ {\e'Y 


K 


= - z + 20 " + + 0 ^ + {^0 




(18) 

(19) 


We ean understand the additional term {^O'^ in Eqs.tfTSl [T9l) as a ‘geometric potential’, 
by analogy to the geometric phase for a spin tha t is ad iabatically transported in a slowly 
varying field (IBerryl. 1 1984 Iwilczek and Shaperd. Il989l) . Since we deal with a second- 


order differential equation, the structure is slightly different from the conventional geometric 
phase: one might also call ‘geometric’ the off-diagonal operator L in Eq. tfTTl) . This operator. 



Figure 3. Illustration of 

the non-adiabatic coupling 9'{z) (see Ea. (fl4l i) for selected energies. Dashed: based on 
Airy and Thomas-Fermi approximations to the condensate density, Eqs.(|4l|5]l. 


involving the derivatives 9' = d6*/dz and 9", is called ‘non-adiabatic coupling’ in the 
following. It peaks roughly where the mean-field potential (j)'^{z) crosses the mode energy 
E, as illustrated in EiglU The Thomas-Eermi approximation 9{z) ~ arctan( 2 :/E) provides a 
simple overview in the dense region, for example (magenta dashed in Eigj3]): 


z > 1 : 


9'{z) ^ 


E 

z^ + E^ 


( 20 ) 


The non-adiabatic couplings are thus confined to the ‘condensate border’ z < E and become 
weak as the energy grows. Conversely, for E —)■ 0, the maximum of 9'{z) shifts into 
the dilute region with a scaling in position (height) roughly proportional to —(logl/E)"/^ 
((log 1/E)^/^), respectively as can be checked from the tunnelling asymptotics of the Airy 
function (dashed gray in Eigj3]). 

In the following, we proceed by adopting first the adiabatic approximation where the 
off-diagonal terms proportional to L are neglected (Sec 12.21) . Non-adiabatic corrections are 
discussed in Sec j2.3[ in particular the role they play for the ‘density mode’ v. 
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Figure 4. (left) Bogoliubov phase mode u in the open channel for different energies, using 
the adiabatic approximation. For the ease of comparison with FiglJ] we have plotted shifted 
potentials E — k'^{z). The wave functions u{z) are multiplied by l/vTO- The bump in the 
potential around z = —1 at low energy is due to the geometric correction Black 

dashed; Coulomb tail of the potential —k'^{z), as given in Eg. (1221) . and corresponding regular 
solution (Ea.(l24lll. 

(right) Closed-channel or density modes v, calculated perturbatively from the adiabatic 
approximation u{z). We have shifted the potentials to — E so that the wave functions 

appear at the energy —E, as expected from Eq.©; they have been multiplied by for better 
visibility. Dashed lines: simple adiabatic elimination v ~ —Lu/rdiz). 


In the adiabatic approximation (subscript ‘ad’), the equation for u can be written in the form 

ad , ’2/ ~ (21) 


+ (E - /C^(^)) Mad = -EWad 


where the potential E — hd{z) is given by Eq.lfTSl) for all z (FigSt/e/f)). At low energies, it 
is similar to the lower (blue) curve in FiglU We call it an ‘open channel’ because we have 
E — k‘^{z) < £' as z —>■ oo so that u is an extended wave right at the continuum threshold. 


with a turning point near 2 ; = —E, as shown in the Figure. Appendix B provides some details 
on the numerical calculation of these solutions. 

Deep in the condensate, we find (black dashed line in FigU/e/f)) 


f^iz) > E : 

- k'^iz) ^ - At + Oiz-\ E^z-^, E^z-^) (22) 

2z Az‘^ 

where the first term recovers the approximation dH). The ‘centripetal term’ ~ 1/z^ arises 
from the first correction beyond the Thomas-Fermi approximation (the coefficient Ci = 1/8 
in Eq.®). The higher-order corrections arise from the next-to-leading order expansion of the 
root E"^ + ^^{z) and from the geometric potential [|6*'(^)]^. The Schrodinger equation for 
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Ua.d{z) therefore matches asymptotically with a modified Coulomb problem: 

dV 


+ Vc{z)f = 0 
1 

2 z Az"^ 


(23) 


dz'^ 

Vciz) = 

an equation that replaces Eq.® obtained above with a simpler argument. The required 
solutions are located just at the dissociation threshold of the Co ulomb potential; they are 
know n analytically and are linear combinations of Bessel functions (lAbramowitz and StegunL 
1972 ) (black dashed in Eig^left)) 

j{z) = y/Wz Jq{EV^) 

y{z) = ^Jt^zVo^Es/^) (24) 

The argument E\/^ of the Bessel functions is familiar from the phase of the WKB solutions 
in Eq. tfTOl) . We have chosen a normalisation such that both Bessel-Coulomb solutions 
have the same amplitude deep in the condensate and their Wronskian is equal to unity 
(Abramowitz and 


mptitude deep 
Steguni 19721) 


TT.r- 1 ■ / ■/ -^y dj 

^[j,y\ =jy -yj 

Deep in the condensate, the Bogoliubov mode can therefore be represented in the form 


(25) 


z ^ 00 : 

fiad(^) ~ COS(5ad “ y{z) sin ()ad) (26) 

where ^ is a normalisation. This formula defines the phase shift (5ad = <^ad(-E') of the 
Bogoliubov mode: the reference case ()ad = 0 corresponds to the Coulomb wave j{z) which 
is regular when extrapolated back to the condensate border (at 2 ; = 0 in the Thomas-Eermi 
approximation). According to the asymptotic series of Jo, Yq, the Bogoliubov mode function 
will match the behaviour deep in the condensate we required in Eq. (fT^ above: 

EV^ S> 1 : 

('2 A1/4 

fiad(^) ~ ^ cos{Es/^ - 7r/4 + (5ad) (27) 

provided we choose the normalisation factor ^ = 1 in Eq. (l26l) (see Appendix A[ Eq. dA. 111) 1. 
We recall that cos6*/2,sin6*/2 —)■ l/\/2 in this limit and that v{z) = 0 in the adiabatic 
approximation. 

Note that the ‘centripetal potential’ —l/(4z^) that arises from the first ‘post-Thomas- 
Eermi’ correction (jfi{z) ^ z — 1/(42;^) in Eq.® is significant in this context. Dropping it 
from Eq. (l23l) . the analytical solutions would involve first-order Bessel functions Ji, Yi which 
are phase-shifted by 7r/2 relative to their zeroth-order counterparts. This could have been 
















Excitations at the border of a condensate 


11 


expected from the long-range character of the centripetal potential, on the one hand. On the 
other, it is interesting to realise that one needs the Jq function in Eq. (l2^ to recover the correct 
behaviour of the Bogoliubov modes at low energies, as required by the U(l) global phase 
symmetry of the mean field theory. We discuss the low-energy limit in more detail in Sec j2.4[ 


2.3. Non-adiabatic coupling and density modes in closed potential 

We now take into account the off-diagonal (coupling) terms in the BdG equations and solve 
Eq.(fT^ for the density mode 

d^fiad 




+ K {z)v^ = -LUi,d{z) 


(28) 


This mode influences significant ly the scattering phase shift 5, as we discuss in a companion 
paper (IDiallo and Henkel 12015h . In addition, it also contributes to the spectrum of density 
fluctuations (dynamical structure factor), as illustrated in Sec l3.3l below. 

The Schrodinger operator on the left-hand side of Eq.(l2^ corresponds to a wedge-shaped 
potential well whose minimum is above zero (EigUrzgfit)). If the right-hand side is neglected, 
we therefore do not have any physically acceptable solution that remains finite for z ±oo. 
Numerically, the inhomogeneous equation can be solved straightforwardly by representing the 
second derivative with a finite difference scheme and solving the corresponding sparse linear 
system. The results are shown in Fig^right) and Fig^left). As expected, the density mode 
is localised in the border region and has an amplitude much smaller than the phase mode. 
The ‘local approximation’ v = —Lu/n?{z) (gray dashed line) captures well its tails, but not 
the reduced amplitude of the oscillatory features (where the second derivative is obviously 
significant). 

Some insight into the inhomogeneous Schrodinger Eq.(l28l) may be gathered by 
considering first the eigenvalue problem 

dfVr, 




-t- K\z)Vn = enVn(z) 


(29) 


for the potential well provided by K^(Z). The eigenmodes provide a convenient basis to 
expand v: 

V(z) = J^bnVn(z) (30) 

n 

The coefficients bn are found by projecting Eq. (l2^ onto Vn, using the natural scalar product 

= J dzVn{z)v{z) (31) 

and choosing the normalisation {vn\vm) = ^nm- We get after two partial integrations 

Gr) — 


(32) 
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One key property here is that the source term in Eq.(l28l). Luad(z), is actually a localised 
function (Figl2Ze/t), thick blue) because the differential operator L involves the derivatives 
d'(z) and d"(z) that tend to zero z ^ E (see Eq.(l20l)l. The matrix elements in Eq.d^ are 
therefore given by convergent integrals. The ‘trapped modes’ Vn are illustrated in Yig^left). 





Figure 5. {left column) Solutions to the eigenvalue problem in the closed channel kHz) (thick 
black), shifted by —E (n = 0 ... 3). We also plot the source term Lu (thick blue) at baseline 

+E. 

The density mode Uad('Z) is shown at baseline —E, solved from the inhomogeneous 
Schrodinger Ea. (l28l l in the closed channel. Thick red = direct numerical inversion of the wave 
operator; overlapping with dashed black = expansion into the lowest twelve trapped modes 
according to Eas. (l^[32l i: dashed gray: local approximation, neglecting the second derivative. 
(right) Squared norm ||u||^ (Ea. (l3Tt ) of the density mode as a function of energy (upper thick 
line). This is compared to the squares (Ea. (l32l ll of its expansion into low-lying trapped 
modes (thin lines, odd modes dashed). At the marked energies, the overlap to the ground mode 
is maximal and minimal, as shown in the left column. 


The closed-channel potential R'^{z) is harmonic only in a narrow range around its minimum. 
Hence, the spectrum is non-equidistant due to the linear asymptotes away from the minimum. 
We derive in [Appendix C| the asymptotics ~ E + [7r(?7, + 


Pseudo-Feshbach resonance. The results for the trapped density mode are summarised in 
Fig^right) where we plot the norm of the density mode Dad (defined as in Eq.d^ below) 
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vs. the energy E. To interpret the oscillating features, we suggest an analogy to the so- 
called Feshbach resonances in atomic and molecular scattering. The physics is essentially 
the same: due to non-adiabaticity, different potential surfaces are coupled. The colliding 
system may thus split up and follow different paths which eventually interfere in the output 
(Stiickelberg oscillations). A particularly strong effect occurs when a localised eigenstate in a 
closed channel becomes degenerate with the incoming wave in an open channel. In ultracold 
collisions, this mechanism operates when a differential Zeeman shift brings coupled spin 
states into resonance; the result is a divergence of the scattering length for specific values of 
the magnetic field (Feshbach resonance). 

In our problem, we also have two potentials, open (—fc^) and closed {hf). But there is 
no tuning parameter available to bring the initial wave into resonance with the closed-channel 
eigenvalues: the denominators e„ in the amplitudes (Eq.(l3^) never cross zero. One can 
even derive the stronger bound e„ > E + vq + ^Jv^ ~ E + 1.22 from the ground state of 
the harmonic approximation to the Hartree-Fock potential contained in hf{z) (see Eq.©). 
There is, however, one possibility for a resonantly enhanced density mode. It is not related 
to a matching of energies, but of wave functions. Indeed, for the energy E 0.32, one 
observes a quite accurate matching in shape and position between Lfiad and the ground state 
^0 (FigStop left)). This leads (as in a Franck-Condon argument) to the strong peak in the 
norm of the density mode 

||i}||^ = (FlF) = Jdzv‘^{z) (33) 

as can be seen in Yig^right) where the probabilities \bff are plotted as a function of energy 
E and compared to the norm ||i)||^. A resonance with the first excited state at ~ 1.5 is 
visible because at a slightly lower energy {Fig^bottom left)), Lfiad becomes orthogonal to 
the ground state vq. At this ‘anti-resonance’, the derivative term in L is significant. 

We have observed that the shape of the closed-channel potential rf‘{z) is relatively stable 
as the energy E increases (compare FigjS^left, top and bottom)). The overlap therefore 
changes chiefly because the turning point and the nodes of the open-channel solution u shift 
with E, as we saw in FigHtZe/r). The other reason is the shifting and broadening of the non- 
adiabatic couplings 6 ', 9" that are involved in the operator L (recall Fig©. 

2.4. Low-energy behaviour 

It is well known that when the Bogoliubov spectrum is continuous, it is gapless and that the 
amplitudes u and v approach the shape of the condensate in the limit E ^ 0. This translates 
the Goldstone mode arising from the global phase invariance (U(l) symmetry) of the Gross- 
Pitaevskii equation @. In this low-energy limit, the phase-density representation of Eqs.® 
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becomes exact. For the sake of simplicity, we stay in the adiabatic basis (fT3]) and get in the 
leading order the BdG equation for the phase mode 

d^-u 


E ^ 0 : 




+ — z)u = 0 


(34) 


This is solved by the condensate 0(z) itself. (For an illustration, see the E = 0.1 curve in 
FigU/e/t).) In the same limit, there is only the trivial solution u = 0 for the density mode. 
We fix the normalisation by continuity with the low-energy limit of the Bessel-Coulomb wave 
that is proportional to the Thomas-Fermi condensate: 

<E^z^l: 

u{z) = ydf(f>{z) ^ j{z) = \A^(1 + 0{E‘^z)) (35) 

We indeed find that the phase shift S{E) is very small in this limit so th at the other Bessel- 
Coulo mb wave y{z) (see Eq.dl^l has negligible weight at low energies (IDiallo and Henkel . 
20151). The spatial range 1 < z -C 1/E^ where this behaviour is relevant opens up wide for 
E ^0. 


3. Applications 

3.1. Equilibrium correlations 

It is well known that the Bogoliubov-de Gennes modes provide a convenient expansion of the 
field operator 

■dE 


/cxh/ 

—= \uE{z)aE + VE{z)a}E} 


(36) 


where the operators {oe) create (annihilate) an elementary excitation with energy E. 
We have assumed a c-number condensate (Bogoliubov shift) for simplicity and added the 
subscript E to the mode functions for clarity. Since the inhomogeneous potential is ‘open’ 
on the dense side, the energy spectrum is continuous. (The integration measure dE/y/ii 


arises from the normalisation of the u, v, see Appendix A ) In thermal equilibrium, we have 
(oe) = 0 and the Bose occupation number 


(aW) = N{E)6{E - E') = 


5{E - E') 


(37) 


qE/T _ I 

because the expansion (l3^ provides a quadratic approximation of the second-quantised field 
Hamiltonian. The elementary excitations contribute even at zero temperature (‘depletion’) 

We focus in the following on low 


because of the operator that appears in Eq. 
temperatures and leave aside the problem of ‘quas i-condensation’ and se l f-consistent mean 
field t h eories in low dirnension s; see for example 
i2002 ): [Mora and Castin ( 2003 ). 


Andersen et al. 


(l2002h : IAI Khawaja et al. 
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3.2. Field correlation spectrum 

Matter-wave interferenee experiments are sensitive to the dynamie field eorrelation funetion 

G{x, y, r) = t + r)) (38) 

where the time dependenee arises from OEit) ~ in the Heisenberg pieture. (Reeall that 
we have set the zero of energy at the ehemieal potential.) Inserting Eq. (l3^ and taking x = y, 
we get the well-known expression 

G{z,t) = |0(z)p 

OO 

+ /— {4(^)^(^) + N{E)) (39) 

0 

We show in Figl^a eontour plot of two terms: the ‘partiele speetrum’ u^iz) and the ‘hole 



-5 0 5 10 15 

position 


Figure 6. Spectral representation of the field correlation function G{x, r). (left column) Focus 
on positive energies: ‘particle mode’ u\(z)', and negative energies: ‘hole mode’ v'^(z). (right) 
Zoom into the low-energy region; note the change in color scale. Dashed and dash-dotted 
white: nodal lines explained in the text. 

spectrum’ Ve{z) (with the energy scale flipped). We recognise in the upper left quadrant 
(particles outside the condensate) a straight nodal line E ^ 2.34 — z (dash-dotted) that is 
characteristic for the Airy function ^^(z) ~ Ai{—E — z). As the modes enter the condensate, 
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his has been regularis ed by introducing 


the nodal lines shift to the pattern £'\/2z 2.405, 5.520 ..., the first few roots of the Bessel 

function Jq (dashed). This approximation is based on the asympt otic form (fT^ and works w ell 
beeause the Bogoliubov phase shift is small: |5(T^) | <^ 7r/2 (see lPiallo and Henkell (l2015h l. 

At low energies (bright central region in Vig^right)), the spectra for both particles and 
holes eonverge to the same limit (see Eg. (1^1 that is essentially given by the eondensate 
density (see See l2.4l) . The oeeupation number N ~ T/E in Eq.d^ thus leads to an infrared 
divergenee of the average densi ty 77 .( 2 :) = (j ( z,0). ’ 
the quasi-eondensate eoneept (iKagan et al.l. I 2 OOOI: 
mainly arises from phase fluetuations whieh ean be subtraeted. See the diseussion of the 
density correlations below. 

Returning to the hole mode ve{z) (Eig^bottom left)), we see that it is eonfined to the 
dense region and follows similar nodal lines as 77 ^; (z) as expeeted from the boundary eondition 
Eq. lfT^ . The eontour plot provides a representation of the so-ealled depletion density 

ro 

■dE 


Andersen et ah . 


200211 : the divergenee 


n,{z) = f—vUz) 
J TT 


(40) 


whieh is simply the zero-temperature limit of the non-condensate density in Eq.(l39l) (a 
measure of quantum fluetuations). We have eheeked that this integral matehes in the dense 
region (i.e., Es/^ ::§> 1) with the eorresponding result for a homogeneous system 0 (loeal 
density approximation with eonstant eondensate f) 


r dk 

J ^ 




■ dk 

2TiEk 


(41) 


where the modes are labelled by the wave veetor k and the dispersion relation is approximately 
linear Ek ~ \/2 kcj). The logarithmie infrared divergenee ean also be eured with suitable 


subtractions (lAndersen et al.L 120021: lAl Khawaja et al.L 120021: [Mora and Castin 


I I 2 OO 3 I 1 . 


3.3. Density fluctuation modes 

As a seeond applieation, we eonsider the leading order Bogoliubov eontribution to the 
dynamic density correlations 

S{x, y, t-t') = \{{p{x, t), p{y, f)}) - n{x)n{y) (42) 

The curly brackets denote a symmetrised operator product for the partiele density p{z, t) = 
fl^z, t)ip{z, t). Its average n{z) = {p{z, t)) does not depend on time (see Eq.(l39l) for r = 0). 

+ A useful parametrisation for the Bogoliubov amplitudes in a homogeneous system is Uk = cosh(?7fe/2), 
Vk = — sinh(pfc/2) with sinhpfc = The dispersion relation is = fc'*-f 2fc^(3^. Therefore in the dense 

limit » Ek: vl Ri e'"“/4 r: (fl/2Ek. 
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The expectation value is worked out using the Bogoliubov shift (13^ and expressed in terms 
of the occupation number s (l37]). using the W ick theorem (gaussian statistics). Our result is 

( 2008 1 and reads 


consistent with Eq.(52) of 


Eckart et al. 


S{x, y, t) = Re {(3^(x, y, -r)A(a;, y, r)} 

OO 

+ (t>{x)(j){y) J— cos{Et) {fE{x)fE{y)N{E) 

0 

+VE{x)fE{y) + {x ^ y)} 

+ 4th order terms (43) 

where the first line involves the correlation function of Eq. (l3^ and the field commutator 

A{x,y,t -f) = \'ijj{x,t), f\y,t') 


— {uE{x)uE{y) e 

TT ^ 




(44) 


Due to the completeness relation of the BdG modes, this goes over into 5{x — y) when f —)■ 0 

|), we use the ‘sum mode function’ 

(45) 


(see [Appendix A| ). In Eq.d 

fE{z) = Ue{z) + Ve{z) 
which is, by a property of the BdG equations, orthogonal to the condensate (f){z) with respect 
to the scalar product (ISTl) . The ‘4th order terms’ of the last line arise from products of four 
Bogoliubov operators oe and Og. Note that the second line features , for a; = y, an i ntegra l 
that is similar to infrared - regularised therm a l dens ities introduced by 


Al Khawaia et al. 

(2002) 

Mora and Castin 

(I 2 OO 31 


Andersen et al. 


( 1200211 : 

This illustrates the consistency of these 
procedures, since their goal is to eliminate from the density spurious contributions attributed 
to phase fluctuations. 

We focus for illustration purposes on the ‘beating’ between the condensate and the 
elementary excitations and show in Eig|7]the local spectrum 


Sheati^Z, E^ 


f^iz) 


TT 


fUz){2N{E) + l) 


(46) 


We find this formula by including the part (l)'^{z)A{z,E) of the first line in Eq.(l431) 
that is proportional to the condensate density. The contour plot shows that the density 
fluctuations are peaking near the condensate border. This is as expected because deep inside 
a (quasi)condensate, such fluctuations are penalised by the self-interaction energy. 

The density fluctuation spectrum does not show any infrared divergence because for 
small E, the sum mode fEiz) behaves proportional to (adiabatic angle 6 ^ 7r/2 in Eq.(fT3])) 

E 

(cos 29 — sin 40) Ue ~ —Mr 
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-5 0 5 10 15 

position 



Figure 7. Spectrum of local density fluctuations due to Bogoliubov excitations in thermal 
equilibrium. We plot Ea.(l46ll which is the essentially the Fourier transform of the second line 
in Eq.(l4^ for temperature T = 1 (using the units specified in Table[TJ. This can be understood 
as a local dynamic structure factor; it is symmetric in energy (only E > 0 is shown), (left) 
Contour plot with nodal lines (dashed) as in Figlh] {right) Cut through the positions z ~ 1 
(peak value of structure factor) and z = 0 (border of Thomas-Fermi condensate). Dashed: 
Bogoliubov amplitudes calculated from ‘phase mode’ u only, ‘density mode’ v omitted. 


for-(lf{z):$>E (47) 

The scaling linear in E at finite temperature can be seen in FigUlright). This plot also 
illustrates that the ‘trapped mode’ ve{z) which is localised near the border (FigJUng/tf)), 
gives a significant contribution (compare dashed and solid lines in FigUlright)). The 
comparison yields the interesting result that the beating between this mode and the condensate 
is actually reducing rather than enhancing low-frequency density fluctuations (set of lower 
curves for z = 0). 


4. Conclusion 


The elementary excitat ions of a Bose condensa t e (Bog oliubov spectrum) are well-known in a 


for a 


harmonica 


1997 


Stringarii 


ly tra pped gas (lAl Khawaja et all 


. 2002: 

Stringari. 

1996; 

Ohbers et al. 


19981). We have analysed in this paper the border region where the 


condensate density smoothly goes to zero, providing a detailed look at the physics beyond the 
homas-Fermi approximation. Previous work has focused on the condensate kinetic energy 


(IDalfovo et al.Lll996l: 


LundhetaL, 


19971), ignoring the contribution of elementary excitation. 


and on the stability with respect to vortex fo r matio n , taking into a ccount motion parallel 


to the border of the condensate dLundh et al.L 


1997 


Anglinl 120011) . The mode functions 
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provided here typically extend into the bulk of the condensate and would correspond in a 
three-dimensional isotropic trap to radially symmetric (angular momentum I = 0) modes. 
Our main result is that the gradient in the condensate density couples elementary excitations 
that are mainly ‘phase-like’ and ‘density-like’, an effect clearly beyond the local-density 
approximation. This leads to density fluctuation modes that are localised near the border 
of the condensate. These fluctuations may be detected with scattering experiments using a 
focused probe bearn that probe the dynamic structure factor locally, similar to the setup of 


Onofrio et al. 


(120001) . Alternatively, one may directly analyse density-density correlations 
when an elongated system is imaged. This may be complemented by launching, with 
a suitable pulse sequence, an elementary excitation coming fr om the bulk (dense quasi¬ 


condensate), similar to the suggestion of 


Brunello et al. 


(120001) . We also believe that the 


methods developed here provide a stepping stone towards a self-consistent description of an 
lomogeneous Bose ga s at finite temperature, using for example the modi fied Popov theory 


m 

of 


Andersen et al. 


(120021) or the Bogoliubov theory for quasi-condensates of iMora and Castin 
(120031) . The border region where the density drops is particularly inter esting here because o f 
the possibility of ent ering a strongly correlate d phase, see for example 


Trebbia et al. 


( 2006 ): 


Jacqmin et al.l (1201 11) . and IVogler et al.l (l2013h . 
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Appendix A. Wronskians and normalisation 

We start by a generalisation of the Wronskian for the BdG problem (two coupled equations). 
Our proposed definition is 

W[u^ V, Ml, Ml] = W[u, Ml] -f W[v, Ml] 

= uu'i — Ulu' + vv[ — viv' (A.l) 

where the prime denotes the first derivatiye. We assume that all modes including the 
condensate f are real. The adyantage of this combination are the following manipulations 
that can be applied to the pair of BdG equations 

Eu = — u" + Hu -I- (ffv (A.2) 

—Ev = — v" + Hv + 0^M (A.3) 

where H = V -|-2|0p — /lis the Hartree-Fock potential. Consider another pair of solutions 
Ml, Ml that solyes the same equations with energy eigenyalue Ei. Multiply Eq. (IA.2l) with 
Ml, and Eq. (IA.3l) with mi, take the sum, and subtract the corresponding equation for mi 
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multiplied by u etc. On the left-hand side, we get {E — Ei){uiu — viv), proportional to 
the integrand of the generalised L^-scalar product in the BdG space. On the right-hand side, 
we find the derivative of the Wronskian W[u,v,ui,vi]: the Hartree potential drops out as 
in the Schrodinger equation; also the coupling terms involving the condensate are cancelled: 

+ viu — uvi — vui) = 0. Integrating and using the boundary conditions (01) for 
z —)■ —cxD, we get 

J dz [ui{z)u{z) - vi{z)v{z)] = Jim ^ ^i] 

In other words: the scalar product can be analyzed locally from the asymptotic behavior of 
the mode functions. The orthogonality between modes with different energies in the discrete 
spectrum follows immediately (the Wronskian vanishes at both ends). 

We continue by analyzing the continuous spectrum for the linear potential. Recall the 
asymptotic form deep in the condensate from Eq. dT^ : 


z -|-cx) : 

u{z) cos{Es/^ — + 5) 


(A.5) 


and similarly for v{z) with amplitude B. In this limit, the two amplitudes are given by the 
rotation back from the adiabatic basis 



Acos9/2 
— ^sin6*/2 


tan 9 = 


<i>\z) 


(A.6) 


where we have used that only the adiabatic mode u ‘survives’ and has amplitude A. (v is 
localised in the Hartree-Fock-like well near the border, see Fig^right).) The scattering phase 
S is therefore the same for both modes u and v. Similar amplitudes Ai and Bi apply for the 
other solution at energy Ei. The Wronskian then becomes (denoting (p = Ey/^ - 7r/4 -f 5 
and similarly for Lpi) 

W[u,v,ui,vi] 


lim 

2^00 


E-Ei 



AAi + BB 


cos cpsmcpi 



(A.7) 


Contributions from the derivatives dA/dz, dB/dz would vanish like Ijz^B relative to this 
term, see Eq. (l20l) . It is natural to interpret this as distribution with respect to the energies 
E, El, to be put under an integral. The trigonometric functions can be re-written into 
sin((^ -f ifi), this gives the expression 
AAi + BBi 


2\/ EEi 


sm 


[{E -f Ei)\/2z — TT/2 -\- 5 -\- ^i] 


(A.8) 
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Since both energies are positive, this is an oseillating funetion as z —>■ cx). It averages to zero if 
integrated over some interval AE > 27r/\/^ and therefore vanishes in the distribution sense. 
The Wronskian is thus given by the phase differenee term ~ sin(<y9 — (fi) 

Wlu,v,ui,vi] 


lim 

2^00 


E — El 
BBi] 


' "'2^/EWi 

= + B^)7i 6{E - El) 


E + El sin[(i? — i?i)\/2i' + <5 — (5i] 


E-El 


where we reeognised in the last fraetion an oseillatory representation of a 5-funetion 

sin(xf) 


lim 

t—>-oo 


X 


= Tr6{x) 


(A.9) 


(A. 10) 


and used the eontinuity in energy of the phase shift and of the amplitudes. The latter sum to 
and in Eq. lfT^ . we ehose the normalisation A = 1. This leads from Eq.d 
to the orthogonality 

■dz 


[— {ui{z)u{z) - Vi{z)v{z)) = 6{E - El) 
J 71 


(A.ll) 


whieh is the main result of this appendix. The symmetry transformation u •H- u and E -f-)- —E 
of the BdG problem (IA.21 IA.3I) gives the additional orthogonality relation 


f^{ui{z)v{z)-Vi{z)u{z)) = 0 
J 71 


(A.12) 


As a eonsequenee, we ean easily eheek that the equations 

dz 


a{E) = J(7 I;{z)ue{z) - 'ij7\z)vE{z)'^ 


= /^ [a{E)uE{z) + a\E)vE{z)^ 


(A. 13) 

translate the standard eommutation relation of the field operator [ijj{x)AKy)] = ^A~y) 
an implementation of the eanonieal eommutation relations for the elementary mode operators 

[a{E), a\E')] = S{E - E') (A.14) 

provided the mode funetions u and v are normalised as in Eq. dA.SI) with A = 1. 
The Bogoliubov shift (eondensate in Eq. dS^ f does not ehange this eonelusion. Eor a 
diseussi qn of the zero mode in the BdG problem and the eorresponding operators, see for 
example IMora and CastinI (l2003h . 
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We use a standard differential equation solver for the Painleve transeendent (Gross-Pitaevskii 
equation ([21)). The solution that eonneets to the Thomas-Fermi profile is aetually numerieally 
unstable, and we mateh it around z ~ 3 with the asymptotie expansion ([5]), keeping typieally 
three terms. 

To solve the BdG equation (fTSl) in the open ehannel (mode u) in the adiabatie 
approximation, a standard forward solver is used: initialise with the tunnelling asymptote (HI) 
and eheek that the potential is linear there. Integrate forward until a position zi ^ 1 and 
mateh to a solution of the modified Coulomb problem (|2^ 

u{z) = a{z)j{z) - I3{z)y{z) (B.l) 

The eoeffieients a, 13 are eonveniently ealeulated with the help of the Wronskians 

a{zi) = W[u,y]{zi), /3{zi) = W[u,j]{zi) (B.2) 

using the normalised Bessel-Coulomb solutions defined in Eq. (l24l) . 

At the position zi, the open potential —k‘^{z) may not yet have reaehed its Coulomb 
asymptote Vc{z) (see Eq.dm)). therefore the eoeffieients a, (3 are still slowly varying. The 
Wronskians (IB.21) satisfy a first-order differential equation that ean be derived using the 
proeedure explained after Eq. (IA.3l) . This yields, for example, 

OO 

= /3{zi) + jdz{-kf{z) - Vc{z))u{z)j{z) (B.3) 

^1 

We ehoose the position zi sueh that the following approximation to the open potential is 
aeeurate enough 


z> Z\ ■. 

- k^{z) - Vc{z) 


E^/{Az) 


(B.4) 


^2 + ziz"^ + ^ 2 ) 1/2 + ^ 2/2 

This term arises from the expansion of the root (0^ -f £^ 2 ^ 1 / 2 . other eontributions (post- 
Thomas-Fermi eorreetion, geometrie potential) are smaller. We find that for z\ ~ 15, the 
relative error is smaller than 10“^ for a wide range of energies. We eompute the integral (IB. 31) 
numerieally with the approximation (IB. II) for u. It eonverges beeause the potential differenee 
seales like l/z'^. One ean avoid the evaluation of oseillatory Bessel funetions for large 
arguments by (i) using their asymptotie form and (ii) shifting the integration eontour into the 
eomplex plane after some point max{zi, 2E} on the real axis. In this way, one is keeping 
elear of the braneh eut of Eq. (IB.4l) at z = ±iE. This proeedure now yields the extrapolated 
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coefficients a, /3. The normalisation factor for the wave funetion u is then 1 /+ /9^)and 
the seattering phase shift follows from tan 5 = f3/a. 

The ealeulation of the density mode v is based on the adiabatie approximation for the 
inhomogeneous Schrodinger equation (fT^ . We represent the differential operator in the 
elosed potential ^^{z) on a grid with a finite differenee seheme. The size of the grid is adapted 
to the support of the souree term Lu. Due to the nonzero minimum of the potential, zero 
energy is not in the speetrum of the differential operator, henee the inhomogeneous equation 
is solved by a straightforward matrix inversion. 

Appendix C. Trapped states 

The elosed potential k^{z) has linear asymptotes on both sides (see its Thomas-Fermi 
approximation in Eqs. (IC.3[ IC.4I) below). Physically allowed eigenmodes therefore join into 
tunnelling solutions and oeeur only for diserete eigenvalues e„ (see Eq. (l29l) . not to be eonfused 
with E whieh remains a eontinuous parameter). 



Figure Cl. Spectrum {cn} of trapped states (colored dots) compared to the Bohr-Sommerfeld 
rule (IC.2b (thick black line). The dashed line gives the analytical approximation (1C. 5b . To 
enhance the difference, the leading term (e — has been subtracted from the action (y- 

axis). 


For the numerical calculation of the trapped states in the elosed potential, we use the 


finite-differenee scheme of the preeeding Appendix B and take a standard sparse eigenvector 
solver. Examples are shown in EigjStZe/f column). A eomparison of the speetrum {cn} with 
the familiar Bohr-Sommerfeld quantisation rule is given in Eig JCll Reeall that this rule is 
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based on the aetion integral 

S{e) =^ +Jdzp{z-,e), 

Piz'N) = ^e-K^{z) (C.l) 

where 2 are the left and right roots of e) (also kno wn as turning p oints). The phase 7r/2 
arises from the Langer eorreetion at both turning points (!MessiahLll995t) . The eigenvalues are 
then approximately given by 


S'(e) = 7r(n + 1), n = 0,l,2,... (C.2) 

The aetion integral, eomputed numerieally, is plotted as thiek lines in Fig JCli and a good 
agreement with the numerieally eomputed eigenvalues is found. For the plot, the eolored 
squares mark the pair (e„,7r(n + 1)). To enhanee the differenee, we have subtraeted the 
leading term (e — from the aetion, see Eq. dC.SI) below. 

The dashed lines in the figure show the Thomas-Fermi approximation to the aetion that 
ean be eomputed analytieally and provides a relatively aeeurate estimate. The elosed potential 
is approximated by 

z <d "zz E — z (C.3) 

z > 0 : R'^ ^ z + V E"^ + z‘^ (C.4) 


These formulas are also useful to estimate the position of the left and right turning points zi^ 2 - 
The aetion integral gives |(e — from the region zi.. .0, and the range 0 ... Z 2 can be 
evaluated with the substitution z = E sinh t. Summing the two, we get 

+ arctanh. (C-5) 

The first two terms give with the Bohr-Sommerfeld rule (1C.21) the eigenvalue spectrum 
e„ ~ E + [7r(n + mentioned after Eq. (l32l) . The scaling law e„ ~ illustrates 

the non-equidistant spectrum in this anharmonic well. Eq. dC.SI) captures relatively well the 
numerically computed action (compare dashed and solid lines in Eig lCll) . except at low 
energies where the Thomas-Eermi approximation fails to reproduce the shape of the potential. 
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